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Abstract 


Many  constrained  optimization  algorithms  use  a  basis  for  the  null  space  of 
the  matrix  of  constraint  gradients.  Recently,  methods  have  been  proposed  that 
enable  this  null  space  basis  to  vary  continuously  as  a  function  of  the  iterates  in  a 
neighborhood  of  the  solution.  This  paper  reports  results  from  topology  showing 
that,  in  general,  there  is  no  continuous  function  that  generates  the  null  space 
basis  of  all  full  rank  rectangular  matrices  of  a  fixed  size.  Thus  constrained 
optimization  algorithms  cannot  assume  an  everywhere  continuous  null  space 
basis.  We  also  give  some  indication  of  where  these  discontinuities  must  occur. 
We  then  propose  an  alternative  implementation  of  a  class  of  constrained  optimi¬ 
zation  algorithms  that  uses  approximations  to  the  reduced  Hessian  of  the 
Lagrangian  but  is  independent  of  the  choice  of  null  space  basis.  This  approach 
obviates  the  need  for  a  continuously  varying  null  space  basis. 


1.  Introduction. 


Many  recent  papers  on  optimization  algorithms  use  a  basis  for  the  null 
space  of  a  given  matrix.  In  these  contexts  the  matrix  often  is  the  transpose 
of  an  nxt  matrix  4(x)  whose  columns  are  the  gradients  of  nonlinear  con¬ 
straints  at  a  point  x.  The  required  basis  for  the  null  space  of  A{x)T  forms  an 
nx(n-f)  matrix  Z  such  that  A{x)T Z  =  0.  When  x  and  thus  A  changes  it  is 
desirable  to  have  the  matrix  Z  change  with  it  in  a  smooth  manner.  Indeed,  in 
analysis  of  optimization  methods,  several  authors  (Coleman  and  Conn  (1904), 
McCormick  (1980),  Nocedal  and  Overton  (1983),  Womersley  and  Fletcher 
(1982))  assume  existence  locally  of  a  continuous  function  Z[x)  that  gives  the 
null  space  basis.  Thus  an  important  question  is: 

Question  1  :  Given  n>ffel,  is  it  possible  to  compute  the  basis  z(j4) 
for  the  null  space  of  a  matrix  A 7  as  a  continuous  function  defined 
for  all  full  rank  values  of  AG.W1*1  ? 

(It  should  be  pointed  out  that  the  null  space  itself  and  the  Orthogonal  projec¬ 
tion  onto  it  are  uniquely  defined;  a  basis  for  the  null  space  is  not  uniquely 
defined  and  so  the  possibility  of  a  continuously  varying  choice  of  this  basis  is 
in  question.) 

This  question  was  addressed  in  a  recent  paper  by  Coleman  and  Sorenson 
(1984).  They  point  out  that  in  several  standard  methods  for  computing  Z 
(usually  as  a  by-product  of  the  QR  factorization  of  A),  Z  is  not  a  continuous 

function  of  A.  However  they  also  show  that  given  any  full  rank  matrix  A,  it  is 
possible  to  compute  Z  as  a  continuous  function  of  A  for  all  A  in  some  neigh- 

A  A 

borhood  of  A  although  this  may  not  be  of  practical  value  in  cases  when  A 
depends  upon  the  solution  of  the  problem  being  solved.  Gill,  Murray, 
Saunders,  and  Wright  (1983)  suggest  a  method  which,  given  a  sequence  \At  j, 


e 


generates  Zk  as  a  function  of  and  the  QR  factorization  of  \-i.  With  this 
scheme,  they  show  that  if  At  ■*  A*.  and  2  I !  ~A*  ||£«,  then  the  sequence 

*=i 

\Zk  j  converges.  This  last  method  appears  to  be  a  practical,  if  perhaps  expen¬ 
sive,  approach.  Note  that  that  with  this  approach,  Z  is  no  longer  a  function 
of  A  and  it  is  possible  that  for  distinct  indices  k  and  l,  A t  ~  A  (or  even 
xk  =  xt)  but  Zfc  *  Zt.  Thus  the  question  of  whether  Z  can  be  defined  as  a  con¬ 
tinuous  function  of  A  remains  to  be  answered. 

In  the  next  section,  we  first  report  that  several  results  from  topology 
show  that  the  answer  to  Question  1  is  "no".  This  answer  also  implies  that  the 
QR  factorization  of  a  matrix  cannot  be  given  by  a  function  that  is  continuous 
everywhere.  We  also  give  a  proof  of  part  of  this  result  that  shows,  in  addition, 
that  the  discontinuities  of  z(A)  can  be  numerous  and  widespread. 

Section  3  discusses  the  dependence  of  some  constrained  optimization 
algorithms  on  the  null  space  basis.  Successive  quadratic  programming  algo¬ 
rithms  that  use  the  reduced  Hessian  of  the  Lagrangian  are  independent  of 
the  choice  of  null  space  basis,  but  similar  algorithms  that  use  a  quasi-Newton 
approximation  to  the  reduced  Hessian  depend  strongly  upon  the  choice  of 
this  basis.  We  propose  a  modification  to  these  methods  which  makes  the 
iterates  independent  of  the  choice  of  the  null  space  basis.  This  approach 
obviates  the  need  for  a  continuously  varying  null  space  basis. 

2.  The  possibility  of  a  continuous  null  space  basis. 

Here  we  consider  a  function  that  gives  a  basis  for  all  or  a  part  of  the  null 
space  of  a  matrix,  and  look  at  the  question  of  what  continuity  properties  it 
can  have.  Let  z  be  a  matrix  function  from  Rnxt  to  /?"**,  where  n  >  f  and 
sfin—  t,  with  the  following  properties. 
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Properties  of  a  partial  null  space  basis  function. 

1.  z(A)ei?***  is  defined  for  all  full  rank  matrices  in  /?*xl. 

2.  z(A)  has  full  rank  s  for  all  such  matrices  A. 

3.  Atz  ( A )  =  0  for  all  A  in  the  domain  of  z . 

If  s  -n— t  we  will  refer  to  Z  as  a  null  space  basis  function. 

The  question  of  whether  such  a  function  z(A)  can  be  continuous  is  an 
important  question  in  topology,  and  was  studied  by  many  mathematicians  in 
this  century.  The  full  answer,  found  by  combining  results  in  Eckmann  (1943), 
Adams  (1962),  and  Whitehead  (1963),  is  stated  below. 

Theorem  2.1  Suppose  z  is  a  function  from  /?***  to  RnX *  satisfying  properties 
1-3,  where  n^2,  lsfcn,  lsssn-f.  A  continuous  function  z  satisfying  these 
properties  exists  in  the  following  cases,  and  no  others. 
n  >  2,  t  =n-l,  s  =  n—t  =  1, 
n  =  7.  t  =  2,  s  =  1, 

n  =  8,  t  =  3,  s  =  1, 

n  =  (2a  +  l)2618c,  f  =  1,  s^2*+8c-l, 
where  aa:0,  6  =  1,  2,  or  3,  cfeO. 

An  immediate  corollary  of  Theorem  2. 1  is  that  the  answer  to  Question  1 
is  "yes"  in  the  cases  when  the  null  space  has  dimension  one  or  when  A  is  4*1 
or  Bxl,  and  "no"  in  all  other  cases.  Thus  in  general,  one  may  not  assume  that 
z(A)  is  continuous.  (When  t  =  n-1,  a  continuous  z{A)  can  be  constructed  by 
a  variant  of  the  method  of  cofactors  for  calculating  the  inverse  of  a  matrix. 
In  the  4x  1  case  the  null  space  basis  is  given  by. 

(Zg  dg  (I4 
— tti  —  a 4  ag 
O4  —a  i  — flg 
— ag  a8  —  Oj 

where  aj,  a2,  a3,  and  a4  are  the  components  of  A.  This  form  is  actually 
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related  to  quaternions;  the  null  space  basis  in  the  8x1  case  has  a  similar 
form  and  is  related  to  Cayley  numbers.) 

Since  the  last  n  -f  columns  in  the  orthogonal  matrix  of  a  QR  factoriza¬ 
tion  of  an  raxf  matrix  A  provides  a  null  space  basis,  it  follows  from  the  above 
theorem  that  the  QR  factorization  of  A  with  n>t  +1  cannot,  in  general,  be 
given  by  an  everywhere  continuous  function  of  A.  However,  the  triangular 
factor  and  the  first  t  columns  of  the  orthogonal  factor  can  be  computed  by  a 
function  which  is  continuous  at  all  full  rank  A.  For  example,  Gram-Schmidt 
orthogonalization  of  the  columns  of  A  has  this  property. 

Theorem  2.2  proves  a  stronger  version  of  Theorem  2.1  in  the  cases  when 
n— f  is  even.  We  show  that  any  function  satisfying  properties  1-3  has  many 
discontinuities  in  this  case.  These  discontinuities  must  occur  within  a 
specified  distance  of  any  given  matrix,  and  they  may  occur  at  quite  well  con¬ 
ditioned  matrices.  The  norms  referred  to  are  the  Euclidean  vector  norm  and 
the  corresponding  induced  matrix  norm. 

Theorem  2.2  Suppose  z  is  a  function  from  Z?"*4  to  /Tnx*  satisfying  properties 
1-3,  with  n-t  even,  and  a  given  matrix  AtR**1  has  smallest  singular  value 
at>0.  Then  for  any  a>0,  z  has  a  discontinuity  at  some  matrix  within  the  set 
\A:  1 1  A~A  1 1  So*  -t-aj,  and  the  discontinuity  occurs  at  a  matrix  with  smallest 
singular  value  equal  to  min{at_i,  a}. 

Proof.  Express  the  singular  value  decomposition  of  A  by 

A  =  UHVT- 

where  U  and  V  are  orthogonal  with  dimension  nxt  and  t xt  respectively  and 
£  is  a  txt  diagonal  matrix  with  diagonal  elements  •  ■  ■  a=ffj>0.  Express 

U  in  terms  of  its  columns  as  Cfsju1  u8  •••***].  Now  consider  the  matrix 

function  4(u)  defined  by 


5 


A(u)  =  [uxuz  •  ■  u*  1  u]  diag{au02,  -.at-ua)  ^ ■  (2.1) 

Let  S  denote  the  n— f  +  1  dimensional  subspace 

jue/?1:  uTux  -  0,  i  =  l . t -lj.  Now  define  the  function  F  on  S  by 

F(u)  =  Z.feD1- 

ll«tf(«)),ll 

where  the  superscript  1  again  denotes  the  first  column.  Note  that  since  z{A) 
has  full  rank  the  denominator  is  never  zero,  so  ||/,'(u)  ||  =1  for  all  ueS. 
Note  also  that,  since  the  columns  of  z(A)  are  orthogonal  to  the  columns  of  A 
for  any  A,  and  thus  to  the  columns  of  U,  F(u)€.S  for  all  ueS. 

Thus  the  function  F  maps  the  unit  sphere  in  S  into  the  unit  sphere  in  5. 
According  to  a  well  known  theorem  in  topology  (see  for  example  Hilton  and 
Wylie  (1962))  any  map  of  the  sphere  in  an  odd  dimensional  space  into  itself 
such  that  F(u )  is  orthogonal  to  u  for  all  u  in  the  sphere  cannot  be  continu¬ 
ous.  Thus,  if  7i -t  is  even,  any  function  z(A)  satisfying  Properties  1-3  will 

have  a  discontinuity  at  A(u)  for  some  choice  of  u  with  ||u  ||  =  1.  Now  note 

that  1 1  A(u)-A  !|^o-t+a  and  that  the  smallest  singular  value  of  A(u)  is 
minjat_i,  aj.  ■ 

This  result  is  fairly  sharp  in  that  it  shows  that  z(A)  must  have  a  discon¬ 
tinuity  in  any  ball  around  A  of  radius  greater  than  <7(,  while  it  is  possible  to 
construct  a  continuous  z{A)  on  any  ball  around  A  with  radius  less  than  at. 
This  fact  follows  from  a  remark  of  Eckmann  (1942).  The  following  theorem 
shows  the  construction  of  such  a  function. 

Theorem  2.3.  Let  AzIF*1  have  smallest  singular  value  >  0,  and  let  Z  be 
any  null  space  basis  for  A.  Then  the  function 

z  (A)  =  [I-A(ATA)-lAT]Z,  (2.2) 

is  a  continuous  null  space  basis  function  over  all  A  such  that  \\A—A\\  <  at. 
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Proof.  Obviously.  ATz(A)  -  0  for  all  A.  To  show  where  z(A)  has  full  rank,  sup¬ 
pose  that  z(A)u  s  0  for  some  A  and  some  nonzero  Since  the  null 

space  of  I—A(ATA)~lAT  is  equal  to  the  range  of  A  then  Zu  ~  Ay  for  some 
ye/?".  This  implies  that  ATAy  =  ATZu  =  0.  But  then 

\\A-A  ||0fe  IK* -ft  II8 

y  y 

_  yT(ATA+ATA)y 

yTy 

*ol 

Thus  z (A)  cannot  be  rank  deficient  if  )jA-A  ||  <aj.  ■ 

Note  that  the  function  z(A)  does  not  give  an  orthonormal  basis  for  the 
null  space  of  A,  but  to  orthonormalize  the  basis  by.  say,  Gram-Schmidt  would 
preserve  the  continuity.  It  is  interesting  to  note  that  the  function  z(A) 
defined  by  equation  (2.2)  actually  gives  the  null  space  basis  for  men 
minimizes  \\Z-Z\\.  (A similar  construction  of  z (A), 

*(A)  = 

which  was  used  by  Goodman  (1982)  in  deriving  a  Newton  method  for  con¬ 
strained  optimization,  can  also  be  shown  to  be  a  continuous  null  space  basis 
over  the  same  region.) 

While  Theorems  2.1  and  2.2  conclusively  answer  Question  1,  a  relevant 
question  for  many  theoretical  results  in  constrained  optimization  is: 

Question  2  :  Given  a  continuous  function  A{ x)  :  DcU"  -*  FT*1 ,  is  it 
possible  to  compute  the  basis  z(x)  for  the  null  space  of  A(x)T  as  a 
continuous  function  of  x  zD  ? 

Obviously  this  is  a  different  question  than  Question  1.  There  are  clearly  some 
functions  A{x)  for  which  a  continuous  z(x)  exists  (e.g.  A(x)  =  a  constant 
matrix  or  a  diagonal  matrix  function  of  z).  Indeed  Wajewski  (1935)  shows 
that  if  D  is  homeomorphic  to  a  sphere,  and  if  A(x)  has  full  rank  for  all  xeD, 
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then  such  a  function  exists.  However  the  general  answer  to  Question  2  also  is 
"no".  For  example,  using  the  same  notation  as  in  the  proof  of  Theorem  2.2, 

for  a  given  A,  consider  a  matrix  function  of  the  form  A{x)  =  A{x)  as  defined 
by  (2.1).  Then  the  function  given  by  the  first  column  of  z(x)  maps  A*  into 
the  unit  sphere  in  1?*.  Now  if,  as  before,  we  restrict  x  to  lie  in  the  unit  sphere 
in  the  the  subspace  S  we  have  a  map  of  the  unit  sphere  in  that  space  into 
itself  which  must  have  a  discontinuity  if  the  dimension  of  the  null  space, 
n—t,  is  odd.  Thus  the  assumption  that  the  function  z(x)  is  continuous  is  not 
justified  in  general. 

3.  The  dependence  of  some  optimization  algorithms  on  the  null  space  basis. 

One  of  the  main  motivations  for  wanting  a  smoothly  varying  null  space 
basis  arises  in  nonlinearly  constrained  optimization  algorithms  that  make 
use  of  an  approximation  to  the  reduced  Hessian  of  the  Lagrangian.  Consider 
the  problem 

minimize  f(x)  (3 

subject  to  c(z)  =  0. 

where  f  is  a  real-valued  function  on  and  c  maps  ft*  to  /?,  where  t  <n. 
Let  j4(x)  denote  the  nxf  matrix  whose  columns  are  the  gradients  of  the  con¬ 
straints.  Given  an  estimate  xk  to  the  solution  of  (3.1),  let  Zk  be  a  basis  for 
the  null  space  of  A(xk)T ,  that  is  an  nx(n  —  t)  full  rank  matrix  such  that 
A(xk)T  Zk  =  0.  Then  the  reduced  Hessian  of  the  Lagrangian  for  problem  (3.1) 
at  X*  relative  to  the  null  space  of  A(xk)T  may  be  expressed  as 
Zk^zL(xk,Xk)Zk.  Note  that  the  reduced  Hessian  is  dependent  on  the  choice  of 
null  space  basis.  We  will  assume  in  this  section  that  Zk  is  chosen  to  have 
orthonormal  columns. 
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Algorithms  for  solving  (3.1)  which  make  use  of  approximations  to  a 
reduced  Hessian  of  the  Lagrangian  have  been  proposed  by  Coleman  and  Conn 
(1984),  by  Womersley  and  Fletcher  (1982)  and  by  Nocedai  and  Overton 
(1983).  These  algorithms  are  either  a  special  case  of,  or  are  very  similar  to 
successive  quadratic  programming,  and  have  the  following  form. 

Algorithm  1. 

Given  xk ,  let 

h±=-ZkMf'ZkrVf{xk) 

v*=-4fc(4fcr4fc)",c(x*)  (3.2) 

or 

vk  ~  ~ (4k4k  )  *c(xk+A|.)  (3.3) 

xk  +  l=Zt+hk+Vk. 

Here,  Uk  is  an  approximation  to  ZkV2L(xk,\k)Zk.  Note  that  A*  is  a  solution  of 
the  homogeneous  equality  constrained  quadratic  program 

minimize  V/ (xk)Th  +)£A TZk Mk Zk Th  (3A) 

subject  to  Ath  =  0. 

Also  note  that  if  the  "vertical  step”,  vk  is  determined  by  (3.2)  as  opposed  to 
(3.3)  this  class  of  algorithms  is  a  special  case  of  successive  quadratic  pro¬ 
gramming. 

We  are  interested  in  the  dependence  of  the  step  (xk  +  l  -  xk)  on  Zk  and  on 
the  continuity  of  our  choice  of  Zk.  An  important  special  case  is  when 

Mk  =  ZkTBkZk 

where  Bk  is  an  nxn  matrix  which  does  not  depend  on  Zk.  For  example,  Bk 
might  equal  V2/, {xk,\k)  or  the  identity  matrix.  In  this  case  A*  is  the  solution 
to 
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minimize  V/  (zk)Th+YfiT Bkh 

hf.fi* 

subject  to  A£h  =  0, 

and  thus,  in  exact  arithmetic,  ft*  and  the  entire  algorithm  is  independent  of 
the  choice  of  null  space  basis  Z 'k . 

However,  when  one  tries  to  generate  M k  oy  quasi-Newton  updates  the 
step  tends  to  depend  strongly  on  the  choice  of  Zk .  Suppose  the  reduced  Hes¬ 
sian  approximation  Mk  is  updated  by  a  method  such  as  DFP  or  BFGS  so  that 
Mk+\  satisfies 

Mk+i  Z*sk  —  ZTyk, 

where  Z  is  either  Zk  or  Zki.\.  Reasonable  choices  for  sk  and  yk  are 

\~xk)  (3.5) 

and 

yk  =  V£(xfc+1.\fc)-7L(xfc+1-A*  ,\fc).  (3.6) 

Other  choices  are  suggested  by  Conn  and  Coleman  (19B4),  by  Womersley  and 

Fletcher  (1982),  and,  comprehensively,  by  Nocedal  and  Overton  (19B3).  Now  if 

Uk  is  a  good  approximation  to  ZpJ2L(zk,\k)Zk,  and  if  Zk+1  is  very  different 

from  Zk,  even  with  a  rank  two  update  we  cannot  expect  Mk+i  to  give  a  good 

approximation  to  Zjt*.iTV2Z,(xk,Xfc)Zk+1. 

Clearly  it  would  be  preferable  for  a  reduced  Hessian  based  algorithm  not 
to  depend  on  the  choice  of  the  null  space  basis  at  all.  Such  an  approach  is 
possible  if,  at  each  step,  we  regard  the  matrix  Mk  as  implicitly  giving  an 
approximation  Bk  to  the  entire  Hessian.  Note  that  if  Uk  is  an  approximation 
to  ZkKPL(xk,\k)Zk  then  ZkMkZkT  is  an  approximation  to  P^L(xk,\k)Pk, 
where 

Pk  =  ZkZkT  =  I-A(zk)[A(zk)TA(zk))-'A(zk)T 
is  the  orthogonal  projector  onto  the  null  space  of  A(zk)T,  and  is  independent 


of  the  choice  of  Zk . 
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Therefore  we  define 

PkBkPk  =  ZkUkZkT 

and  extend  our  definition  of  Bk  to  approximate  the  rest  of  the  Hessian.  The 
idea  of  a  reduced  Hessian  method  is  to  ignore  the  contribution  of 
Pk^iL(xk,\k)(I-Pk)  to  the  Hessian  so  we  will  take 

PkBk(l-Pk)  =  ( I-Pk)BkPk  =  0. 

Now  the  term  ( I-Pk)Bk(I-Pk )  has  no  effect  on  the  current  step,  but  it  could 
affect  all  parts  of  Bki.  t  and  thus  subsequent  steps.  Just  as  a  scaled  identity  is 
often  used  as  an  initial  Hessian  approximation,  we  will  let 

(I-Pk)Bk(I-Pk)  =  (3k(I-Pk)Hl-Pk)  =  Pk(I~Pk). 

The  scaling  factor  pk  is  an  approximation  to  1 1  VzL(xk,Xk)  |  | ;  for  example 

pic  =  \\Mk  ||  is  reasonable  and  minimizes  the  condition  number  of  1 1  Bk  |  | . 

Thus,  given  Mk,  our  implicit  Hessian  approximation  is 

Bk  =  Pk  Bk  Pk  +Pk  Bk  (I-Pk ) + (/ -Pk  )Bk  Pk  +(/  -Pk  )Bk(I-Pk) 

=  ZkUkZk7+pk(I-ZkZkT). 

Now  given  the  approximation  Bk  to  the  entire  Hessian,  the  reduced  Hes¬ 
sian  approximation  corresponding  to  a  new  basis  Zk  +  l  is  given  by 

Zk+xTBkZk. ,  =  TiMkTk+pk(I-TlTk) 

where 

Tk=Zk7Zk+l. 

This  matrix  may  then  be  updated  to  give  Mk+ 1- 

The  above  logic  leads  to  the  following  algorithm. 

Algorithm  2 

1.  Given  xk  and  Mk  compute  li*  and  vk  as  in  Algorithm  1. 

Z  Let  xk+i  =x k+ht+v/t. 

3.  Compute  Zk+ j,  Tk  =  ZkTZk+i  and  pk. 


4.  Let 
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a  =  n{Mk-pki)Tk^ki. 

6.  Update  Mk  to  get  Mk+{  using  the  DFP  or  BFGS  update. 

with  secant  equation  Mk¥X Zk+\Tsk  =  Zk+Jyk. 

6.  Set  k  to  k  + 1  and  go  to  step  1. 

In  step  5  a  choice  of  s*  and  yk  such  as  that  given  in  (3.5)  and  (3.8)  is  rea¬ 
sonable  as  are  any  of  the  Zk+ 1  based  choices  given  in  Nocedal  and  Overton 
(1983,  p.28).  However  if.  instead,  the  update  involves  ZkTsk  and  Zk Tyk .  as 
does  the  one  used  by  used  by  Coleman  and  Conn  (1984,  p.  780),  then  the 
update  in  step  5  should  come  before  step  3. 

It  should  be  noted  that  this  algorithm  involves  three  extra  matrix  multi¬ 
plications  over  Algorithm  1  for  a  total  of  2(n-f  )3+n(n-f  )2  extra  multiplica¬ 
tions.  If  n-t  is  small  relative  to  n  this  is  much  less  than  the  work  involved  in 
the  computation  of  Z;  it  is  also  less  work  than  the  modification  proposed  by 
Gill.  Murray.  Saunders,  and  Wright  to  make  Z  vary  smoothly,  which  was  men¬ 
tioned  in  Section  1.  The  benefit  of  this  extra  effort  is  that  the  algorithm  is 
independent  of  the  choice  of  basis  as  described  below. 

Theorem  3.1.  If  Algorithm  2  is  used  with  the  initial  matrix  Mq  =  ZqTBoZo 
where  B0  is  independent  of  Z0,  and  if  sk  and  yk  are  given  by  (3.5)and  (3.6), 
then  the  iterates  produced  are  independent  of  the  choice  of  Zk  for  k&l. 

Proof.  One  need  only  note  that  Algorithm  2  is  equivalent  to  the  following  algo¬ 
rithm  if  we  let  Bk  -  Zk  Mk  Zk  r. 

Algorithm  3 

1.  Given  xk  and  Bk  compute  hk  by 

minimize  V/  (xk)Th+}ffiT  Bkh 
heir' 

subject  to  A%h  =  0. 

Compute  vk  as  in  Algorithm  1. 
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2.  Let  xk+l  =  xk  +A*  +vk . 

3.  Compute  Pk+\.  and  (3k. 

4.  Let 

=  Pk+l(Bk—PkPk)Pk  +  l+Pkf- 

5.  Update  Bk  using  the  DFP  or  BFGS  update, 

with  secant  equation  Bk¥i  Pk^iSk  =  Pk+ xyk. 

6.  Set  k  to  Jfc  + 1  and  go  to  1. 

Note  that  this  algorithm  does  not  involve  Zk  at  all;  thus  the  equivalent  Algo¬ 
rithm  2  is  independent  of  Zk.  ■ 

Algorithm  2  is  very  similar  to  the  algorithms  proposed  by  Coleman  and 
Conn  (1984),  by  Womersley  and  Fletcher  (1982),  and  by  Nocedal  and  Overton 
(1983).  We  have  of  course  only  described  it  very  generally,  but  the  main 
difference  is  the  subspace  shift  in  Step  4.  Because  of  this  it  cannot  be  con¬ 
sidered  a  special  case  of  these  algorithms,  and  questions  of  convergence 
need  to  be  considered.  We  believe  that  the  subspace  shift  is  a  good  alterna¬ 
tive  to  trying  to  ensure  some  kind  of  continuous  change  in  the  null  space 
basis. 


l 
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